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SUMMARY 

Computations in tribology and material science at the atomic level present considerable difficulties. 
In this paper computational techniques ranging from first-principles to semi-empirical and their limita- 
tions are discussed. Example calculations of metallic surface energies using semi-empirical techniques are 
presented. Finally, application of the methods to calculation of adhesion and friction are presented. 


INTRODUCTION 

Computational material science or the theory of materials has recently come of age. Calculation of 
properties of real materials at the atomic level such as grain boundary or dislocation energies or the dy- 
namics thereof, which in the recent past have seemed intractable, now have some hope for realistic 
modeling. An even more startling assertion is that modeling of tribological phenomena is also now feasi- 
able. Tribology is a particularly difficult field for theoretical studies at the atomic level because of the 
different possible materials in contact under high loads often with high degrees of disorder, but as 
Landman’s et al. have shown (paper in this proceedings) a great deal of progress has been made in 
approaching problems of practical interest. This paper is concerned with developments which have 
enabled this progress. An effort will be made to define and explain the various approaches for calculating 
the configurational energy and to give an understanding of such approaches. Finally, we will concentrate 
on two methods for determining energetics, the embedded atom method (EAM) and equivalent crystal 
theory (ECT). We concentrate on EAM because it is the method of choice at present and has success- 
fully been used to treat many problems concerning defect energetics. We include ECT, because it is a 
new method which has the capability of treating a wider class of materials with good quantitative agree- 
ment for defect energetics in situations where there is a large deviation from equilibrium. An example of 
how one might calculate surface energy, adhesive energy, and energy as a function of position for sliding 
one metal over another will be given with the intent of encouraging the unfamiliar reader to apply these 
relatively simple techniques to problems of interest. 


BACKGROUND 

First, we address the issue of what has changed in atomistic modeling in order to enable examina- 
tion of problems of interest in tribology. First principles (ab-initio) approaches have dealt with perfect 
single crystals in which lattice periodicity (ref. 1) reduces the calculation to a single unit cell. The results 



of these calculations were impressively successful, in that accurate agreement with experiment was 
obtained in predicting the band structure, transport and magnetic properties of solids. Properties of 
interest to the material scientist, however, are generally related to “defects” in atomic structure which 
involve a partial loss of periodicity and thus require calculations over a large number of atoms, for 
example, dislocations, grain boundaries, and interfaces between different materials. The breakdown of 
periodicity and the involvement of many atoms clearly complicates such problems, for example, we show 
a twist boundary for an fee (100) interface (fig. 1). One can see that the geometry is quite complex and 
that there are many nonequivalent atoms which must be considered in the calculation. Until recently, 
the lack of adequate affordable computing power combined with lack of sufficiently efficient methods for 
performing quantum mechanical (ab-initio) calculations made such complex problems difficult. For ex- 
ample, it only recently (refs. 2 and 3) became possible to predict which simple structure, fee, bcc, or hep 
for an elemental metal had the lowest energy. These limitations lead to the use of alternative approaches 
such as pair-potentials (two-body forces) to calculate properties of defects. An additional problem in any 
approach is that the energies of interest are obtained from the difference in total energy between an ion in 
the solid and in a reference structure. Each of these quantities is a large number whereas the difference is 
a small number thus causing a significance problem in obtaining structural energies. We show a simple 
example (fig. 2) using the jellium model where the charge density of the positive ions in a metal is 
smeared out uniformly over the solid. The energy difference in creating two surfaces arises from changes 
in the surface region. The bulk remains the same, thus, integrating over volume with and without the 
surfaces and subtracting gives a difference in energy small compared to the integral over the entire 
volume. This situation applies regardless of the approach used to calculate defect energies. 

For the next issue, we again refer to figure 1. The rigid twist boundary defect shown is not the 
minimum energy configuration, consequently a complex search must be performed to find the minimum 
energy configuration. Finally, if temperature and dynamic effects are to be included the situation be- 
comes even more complex. Having outlined the problem we now examine some approaches. 


APPROACHES FOR DESCRIBING THE INTERACTIONS 

In this section we will give a general discussion of approaches for calculating defect energies. The 
order will be inverted in that we will start with the most general, proceed to the least general, and build 
to better approximations from the least general. We address four approaches, first-principles calcula- 
tions, pair-potentials, many-body potentials, and semi-empirical methods. The purpose of these dis- 
cussions will be to give a qualitative understanding of the methods and to provide a starting point for a 
literature search should a more comprehensive understanding be desired. 


First-Principles Calculations 

This is the most fundamental approach for calculating material properties. In principle, the only 
inputs are fundamental constants and the atomic number of the atoms of interest and it is not materials 
limited. Recently, techniques have been developed to treat some of the problems discussed, i.e. , loss of 
periodicity, energy minimization, and dynamic effects. The loss of periodicity is handled by a trick called 
“super cells” (ref. 4). A quasi-periodicity is produced in which a large cell which mimics the defect is 
repeated throughout the solid. The hope is that if the cell can be made sufficiently large, the energetics 
will converge to that of the real defect. Of course in highly disordered systems such as one might find in 
tribology this approach will not work as well. Carr and Parrinello (ref. 5) have developed a new tech- 
nique which allows minimization of the total energy as a function of the electronic degrees of freedom and 
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also allows the study of dynamic and temperature-dependent effects. However, first-principles calcula- 
tions remain cpu time intensive. 


First-Principles (Ab-initio) Methods 

We now proceed to give a more detailed description of first-principles calculations. These normally 
involve solving the Kohn-Sham equations within the local density approximation (LDA) (ref. 6). These 
are a set of equations using the mean-field approximation which are written as a set of one-electron equa- 
tions where the electron is moving in the mean field of all of the other particles. The total energy is 
expanded in terms of the electron density, since it can be shown that the total energy is a functional of 
the electron density. In the LDA only terms that depend on the density at a given position are retained 
and not terms in the gradient or higher order derivatives. These equations can be written as 

((-1/2) V 2 ♦ V(r))* k (r) > « k * k (r) (1) 

where and e ^ are the one-electron wave functions and energies and k is the wave number and 

V(r) = *(r) + V xc [/j(r)] ( 2 ) 

where $(r) is the electrostatic energy V xc [p(r)] is the exchange and correlation energy, and p(r) is the 
electron density. Square brackets indicate a functional of the electron density. The wave functions are 
often expressed in terms of some appropriate basis set (ref. 1). The electrostatic energy is given by 

$ (r) = J dr' p(r)/|r - r'| Zj/|r - Rj| . (3) 

The first term is the electron-electron interaction and the second is the electron-ion interaction, r is the 
electron coordinate, R is the ion coordinate, and Z is the atomic number. The electron density is given 
by 


?( r ) " Eocc. M 2 ^ 

where occ. refers to the sum over occupied states. With the geometry and an appropriate approximation 
for the wave functions, these equations are solved self-consistently, that is, one iterates until the output 
and input density or potential agree to within some prescribed limit. The total energy is then obtained 
within the LDA as 


E [ p) - E KeM + E EsM + E Xc(H + E ion- ion ( 5 ) 

where KE refers to the kinetic energy, ES to the electrostatic energy, XC to the exchange and correla- 
tion energy, ion-ion refers to the ion-ion interaction energy, and the brackets indicate that the total 
energy is a functional of the electron number density (atomic units are used throughout). Although the 
first-principles methods are based on the firmest theoretical foundations, there are still many approxi-ma- 
tions used to get a solution, such as the mean-field approximations and the LDA. Of the methods to be 
discussed, only ab-initio methods can be used to calculate the electronic structure of the solids and are 
material independent. 
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Pair-Potentials 


In this section we will start with a more general discussion following Carlsson and Ashcroft (ref. 7), 
which will lead naturally to the following sections. What is the functional form of the cohesive energy 
E coh(^l v • *>Rn) in a system consisting of N atoms with nuclei at sites { R- } ? A simple guess for express- 
ing the cohesive energy then would be a cluster expansion in pair and higher-order interactions 


E coh = E 1 + E 2 + 


= E, 


E 3 + . . . 

+ (l/2) £ ' V 2 (R i E j ) + (1/6) £ ' VgfRiRjRj + . . 

iJ hiM 


where the prime indicates that no two indeces have the same value and represents external forces. 
The pair-potentials refers to keeping only the V 2 term in the expansion which dominates when the 
higher order terms are small in magnitude compared to the pair term. In practice, however, the appli- 
cation of pair-potentials have not followed from such an approach, but by simply assuming some form for 
V 2 and then empirically fitting to some physical parameters, for example, a Morse potential 



The unknown parameters could be selected to give the correct cohesive energy, bulk modulus and equi- 
librium lattice constant, etc., depending on the number of fitting parameters in the expressions. Once 
these parameters and the defect geometry are established, the energy is calculated by performing a sum 
over all pairs of atoms. For some solids such as rare gas solids this approach might be sufficient. In 
metals, where volume-dependent terms are important (ref. 8), it is not. For example, in figure 3 we show 
the results of a first-principles calculation for a grain boundary using the jellium model (ref. 9) where the 
grain boundary is represented by a lower jellium density. We can see that the electron gas has rear- 
ranged from the uniform gas. This term would be omitted in the pair-potential approach. There are 
other deficiencies. For example, they require that the elastic constants satisfy C 12 = C 44 which seldom 
occurs in real metals. In addition, they require that the vacancy formation energy equals the cohesive 
energy which is also not the case. The vacancy formation energy is typically of the order of one-third of 
the cohesive energy (refs. 7 and 8) for metals. The principal advantage in pair-potentials is simplicity, 
however they may be useful in predicting trends when stearic effects dominate. 


Many-Body Potentials 

For completeness we now give an example of including the term in the series given in equa- 
tion (6). Stillinger and Weber (ref. 10) have proposed a potential energy function for Si which is con- 
structed to account for the strong directional bonding present in such elements, since pair-potentials had 
failed to predict behavior for large excursions from equilibrium. They propose for the pair term 
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V a W = 


a(t p - r“ q )exp((r - a)" 1 ), r 


< a 


( 8 ) 


0, r > a 


and for V, 


where 


v 3( r i» r j’ r k) ^( r ij’ r ik’^jik) + ^( r ji’ r jk’^ijk) + * 1 ( r ki’ r kj’^ikj) 
h ( r ij> r ik^jik) = Acx p(t ( r y - a)' 1 + 7 (r ik - a)' 1 ) (cos 0 jik + 1/3) 


(9) 


( 10 ) 


where h vanishes identically for r > a. The “ideal” tetrahedral angle 0 t is such that cos 0 fc = (-1/3) 
so that the trigonometric part of the expression discriminates in favor of pairs of bonds emanating from 
vertex i with the desired geometry. The fitting parameters were selected through a search which insures 
that the diamond structure is indeed the most stable arrangement of particles, and that the melting point 
and the liquid structure inferred be in reasonable agreement with experiment. For a further discussion of 
the inclusion of higher order terms see references 7 and 8. 


Semi-Empirical Methods 

This category of approaches tackles the many-body problem in a different fashion. The term semi- 
empirical refers to determining a functional form for the cohesive energy based on some physical model, 
but including some parameters determined by fitting to experimental properties. Once these constants 
are determined, the form is used to calculate the energy or dynamic behavior or other properties of in- 
terest such as defect energies. This approach has proven to be extremely effective in accomplishing its 
goals. This section concentrates on three methods, the method of Finis and Sinclair (refs. 8 and 11) 

(FS), the embedded atom method (refs. 12 to 14) (EAM), and equivalent crystal theory (refs. 15 to 17) 
(ECT). The next section concentrates on giving examples from EAM and ECT. The reasons for this 
idiosyncratic choice is that the authors have no experience in using (FS), however have used EAM and 
are co-developers of ECT. ECT treats a wide class of materials with a high degree of agreement with 
experiment and first-principles calculations. EAM has been applied to a wide variety of problems with 
considerable success. 

FS. - Finis and Sinclair (FS) express the total energy of an assembly of atoms in two terms 

U,„, - U N * U p (11) 

where is the N-body “cohesive” term and Up is a conventional central pair-potential repulsive 
term. is selected to have the form 
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( 12 ) 


and 


where 


U N A E f (pi) 

i 

Up - (1/2) £ v(r m ) 

i, j 

P, - E*K) 

j 

R ij = Kl = K - R j| 


(13) 


(14) 

(15) 


The ^(Ry) are the electron charge density at site i from an ion atom located at site j and V(R ■) is a 
pair repulsion between atoms at site i and j. Therefore, equation (11) is interpreted as the sum of a 
volume dependent term, U N , and a pair wise summation U p . The term f(/o) is chosen to be \t p in 
order to mimic the result of tight-binding theory (ref. l). The unknown constants are obtained from fit- 
ting to the cohesive energy, elastic constants, and the bulk modulus. The reader is referred to reference 8 
for implementation of the model. Application of the model is similar to methods to be discussed next. 

Universality of binding energy relations . - We digress for the moment to discuss an issue which will 
be important in the discussions of EAM and ECT. Rose, Smith, and Ferrante (refs. 18 and 19) (RSF) 
discovered that there were similarities between certain binding energy relations, i.e., binding energy 
versus distance. They found that the functional form of the scaled binding energy curves describing 
cohesion, adhesion, chemisorption, and diatomics for the case of no charge transfer was the same to a 
high degree of accuracy (fig. 4). Stated more formally, we write the energy versus separation as 

E(r) = AE E*(aj ( 16 ^ 

where 

- (R - R,)/< (”) 

AE is the binding energy, a is a scaled length, R is the distance between particles, R e is the equili- 
brium distance and f is a scaling length. RSF found that a simple functional form was an accurate 
representation of the function E (a ) first proposed by Rydberg for diatomics 

E*(a*) = -(l + a*) exp ( -a*) ( 18 > 


(19) 


This result was dubbed the “universal binding energy relation” (UBER). The application of this result to 
EAM and ECT will be discussed in the next sections. 


and the length scaling was chosen as 


= (AE/(d 2 E/dR 2 ) R J V2 
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EAM * -The embedded atom method (EAJvl) can be considered as a breakthrough in atomistic 
simulations of metals. It has been successfully applied to problems such as phonons, liquid metals, 
defects, fracture, surface structure, surface segregation, etc. (ref. 20). Although it is very similar in 
structure to FS, the justification for it is different. It is based on the effective medium theory of Stott 
and Zaremba (ref. 21) or the quasi-atom approach of Norskov and Lang (ref. 22) in which the energy to 
embed an atom in jellium is calculated. The total energy of the system is written sis 

N 

E.oi - E E i ( 2 °) 

i 

where 

E i = F iki) + (1/2) £ *(R U ) (21) 

j * i 


and /> h>j is the total electron density at atomic site i due to all the other atoms in the system, F- is 
the embedding energy for placing an atom at site i and $ ( R..) is a short-range pair-interaction rep- 
resenting the core-core i and j separated by R^. N is the total number of atoms 

in the system. The host electron density at site i is approximated by the linear superposition of all of 
the other atomic densities 

- e «;K) (22) 

j 


where p j |Ry j is the atomic density of atom j at the distance R- from its nucleus and 

?j( R ) = *W R ) + n W R ) ( 23 ) 

where n g + n d , the total number of s and d electrons, is constrained to be 10 for Ni, Pd, and Pt; and 
it is 11 for Cu, Pd, and Au. The core-core repulsion term takes the form 

*K) ■ z iK) z XK ( 24 > 

where Foiles, Baskes, and Daw have parameterized a simple function for the effective charge Z(R) as 

Z(R) = Z 0 (l + /3R V ) exp(-aR) ( 25 ) 

or ^er to guarantee the short range for the pair interaction, and Zq is the nuclear charge. Originally, 
Foiles, Daw, and Baskes (FDB) used a spline fit to define the embedding energy F(/>), but later found 
that defect energetics could be better approximated by requiring that for the case of the cohesive energy 
with an isotropic expansion or contraction of the solid that F(/>) be required to agree with the UBER 
(eq. (18)) giving 
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F ( R ) = Ecohi 1 + wt)exp(-wc) - <j> 


(26) 


where $ is the second term in equation (21), E coh is the cohesive energy, x = (r/r Q - 1) and 
- = 3 (V 0 B 0 /E coh )‘/ 2 . V„, B q , and r Q are the equilibrium volume, bulk modulus, and Wigner-Seitz 
radius, respectively. The unknown constants are again obtained by fitting to such properties as the 
vacancy formation energy, elastic constants, and diatomic parameters. It is unnecessary to perform the 
fitting for many materials since they are provided by FDB in their papers. Another point to emphasize is 
that once the embedding energy is specified for a given metal the same function is then used for any 
defect. Alloys have been treated in EAM (ref. 23) by keeping the same embedding function for a given 
metal and constructing the density by the overlap of atomic densities for all species present. For exam- 
ple, in a Ni-Cu alloy, the geometry would be specified, then taking a Ni site, the overlap contribution to 
the electronic charge density from both species at the Ni site is calculated. Finally, that density is used 
in the Ni embedding energy function to determine the contribution for the Ni atom at that site. A 
geometric mean of the two pair repulsion terms is used for the interaction between different atoms, i.e., 

*ab( r ) = [*aa( r )*bb( r )] 1/2 - W 

Johnson (ref. 24) has proposed an alternate form for the pair-interaction which gives the same invariance 
with respect to electron density transformations as monatomic models. 

Equivalent crystal theory (ECT) . - Equivalent crystal theory is a new technique developed by 
Smith, Banerjea, and co-workers (ref. 25) which also uses the UBER in order to implement the method. 
The technique gives quantitatively accurate results for the surface energy and surface relaxation of metals 
and covalently bonded solids. The application of EAM to bcc metals (ref. 26) and to covalently bonded 
solids is limited (ref. 27) and thus ECT has fewer material limitations than EAM. The idea is based on 
the fact that the UBER gives the ground state energy relation for a given solid. Consider first a single 
crystal of an elemental solid. Next introduce a defect or an array of lattice defects into the crystal. The 
total energy E(defect) of the crystal containing the defects is equal to energy of a single crystal E(crystal) 
plus a perturbation series, where the perturbing potential is the difference between the array of ion core 
potentials of the crystal containing defects and those of the single crystal: 

E (defect) = E (crystal) + Perturbation series (28) 

E (crystal) is given by equation (16) where a = (r wg - r wsE)/^ r WSE * s equilibrium Wigner-Seitz 
radius [3/(47rr wgE ) = bulk atom density], 1 = [AE/(127rBr wgE )] 1 ^ and B is the bulk modulus of the 
ground state single crystal. At this point we present the idea of the “equivalent crystal.” The 
equivalent crystal is a conceptually ideal crystal of the material which is either expanded or con-tracted 
from the ground state crystal and thus its energy is given by the UBER at some different scaled distance, 
a Q . The procedure is to find the value of a* for which the perturbation series = 0. If this can be 
accomplished then 


E (defect) = AE E*(a,J ) (29) 

which can be evaluated simply by using the UBER. In order to implement ECT, we must have a method 
to represent the perturbation series. As before, the purpose of this presentation is to show how to use 
ECT, hence no formal justification will be given and the reader is referred to references 25 and 17 for this 
aspect. The perturbation series has matrix elements which are integrals over products of the density and 
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the difference in potential between the perturbed and unperturbed state. These integrals are approxi- 
mated by 


g a R p exp(-aR) 


( 30 ) 


for nearest neighbors and 

g « R p exp(-aR) ex p(-R/A) (31) 

for next nearest neighbors and beyond, where p = 2n - 2, n = principal quantum number and A = elec- 
tronic screening length. It is assumed that the form of the electron density is that of the highest partially 
occupied s orbital. Screening is introduced for next nearest neighbors by addition of the term 
exp(-R/A). The calculation proceeds for solving the perturbation equation on a site by site basis, so that 
the energy change is written as the sum over the changes at each nonequivalent site given by 

S 

[Ej(<fefecfc) - Ej(crys£a/)] (32) 

1=1 


where the Ej’s are the energies associated with site 1 for the defect crystal and the ground state crystal, 
respectively. In many defects where a high degree of symmetry is maintained such as a surface, the value 
of s in equation (22) can be quite small. We now present an older version of ECT which is sufficient to 
calculate surface energy and surface relaxation for our example calculation. The more general method 
will be outlined after this discussion. We now introduce the actual working equations for application of 
the perturbation equation (32) for 6E = 0 for a given site I, to next nearest neighbors 


b, R 


1 ( 1 ) 


-aRi(l) 


+ b 2 R£e 


p -aR^l) -R2(l)/A 


E 

defect n 



4 - 


defect nn 


(33) 


where bj and b 2 are the number of nearest (next) neighbors, respectively; and Rj and R 2 are the 
nearest (next) neighbor distances, respectively. We solve this equation for Rj, the equivalent crystal 
nearest neighbor distance, which has a simple geometric relationship to the Wigner-Seitz radius and all of 
the neighbor distances since we are dealing with a perfect crystal. To restate, an atom in a defect in the 
material experiences an environment as if it were in a perfect crystal, which is expanded or contracted 
from its equilibrium state to a new nearest neighbor distance given by Rj. The expression for the sur- 
face energy, for which we will give an example later, is given by 


a = (ae/a) £ 
1=1 



where 
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F (a ) = 1 + E*(a*) = 1 - (l + a*)e a 


(35) 


where A is the surface area, s = number of layers with energy different from the bulk (* 6), 

M = number of nearest neighbors for layer 1, 0 ml = 1 if a ml < 0 and 0 ml = 0 otherwise, 

= ^i(l)/ c i ” r WSE ( c i rat ^° between the nearest neighbor distance and the Wigner-Seitz 

radius), a^j = Rmi/cj - r wSE’^ml * s ^he distance between atoms m and 1, and is the number 

of nearest neighbors of atom m or 1, whichever is smaller. The first term in equation (34) gives the 
deviation from the minimum in the binding energy curve and thus the increase in energy from being in 
the defect. The second term is included to deal with the fact that when relaxation is included bonds are 
distorted inhomogeneously, as opposed to the first term which depends on the average nearest neighbor 
distance. The second term is zero when bonds are not compressed. For the rigid surface energy (i.e., no 
configurational relaxation) the second term can be neglected. ECT in this form has essentially only one 
fitting parameter, o, which is obtained from the vacancy formation energy. The other physical param- 
eters such as the cohesive energy, bulk modulus, and equilibrium lattice parameter appear explicitly in 
the model. 

Since the purpose of this presentation is to give an overview of the methods, we will not go into 
details for the modification of ECT so that a wider class of materials and situations can be treated. We 
will only indicate the general procedures leaving the details to reading reference 25. The defect formation 
energy is now written in terms of four components: 











ej = AE " 

F * 


+ f1 



a 3 (i Jjk) 

* E F ' 

p,q) ‘ 






j. k 


p> q 

i. j 


(36) 


This expression is now in the spirit of the many-body expansion presented in equation (6). Note that 
for the case of surface relaxation that equation (34), which is a simpler version of equation (36), suffices. 
In this approach each term is treated as though it is an independent perturbation. The first term is that 
presented in the original formulation of ECT, the second, the bond-compression term, although formu- 
lated differently, has the same purpose as already presented. The third term is included to handle the 
situation where bond angle changes are important as in semi-conductors. The fourth term is included in 
order to give the correct values for the shear elastic constants. In this formulation there are four fitting 
parameters compared to one, a. These four parameters are obtained from the vacancy formation energy 
and the elastic constants. The other input parameters remain as before. The procedure is more com- 
plicated but not difficult, because each perturbation is treated as independent. There are now four 
equations similar to equation (33) but no more difficult than it to solve for the a 5 s. The results for 
EAM and ECT will be shown in a later section. Finally, alloys (ref. 28) are treated by an extension of 
the basic ECT ideas in which the a’s reflect whether the neighbors are of type A or B. 


EXAMPLES OF SURFACE ENERGY CALCULATIONS 

We now give an example of a surface energy calculation for a rigid fee (100) metal surface using EAM 
and ECT to demonstrate the relative simplicity of such calculations. Each atom in each plane is equiva- 
lent to all other atoms in the plane, therefore it is only necessary to solve the EAM or ECT equations for 
a single atom in each plane and the only remaining question is how many planes are needed for conver- 
gence. In EAM, the first step is to pick an atom in the plane, say atom A, and then evaluate the electron 
density at A from all of the other atoms in the crystal (eqs. (22) and (23)). The density is evaluated 
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from analytic expressions (ref. 12). Once the density is known, one determines the embedding function, 
F(p) (eq. (20)), for the given metal at that density. Then one need only evaluate the pair term in 
equation (20) using the parameters given in reference 12 and then evaluate the contribution to the energy 
for this atom from equation (20). The calculation is finished for the top plane. The same calculation is 
repeated for a single atom in the second plane. This procedure is continued in each plane until there is 
no change in the energy from the bulk value. For EAM one would subtract the bulk energy for each 
atom, sum and divide by the area of a primitive cell in order to calculate the surface energy, i.e., 

» - E (e? 1 ”' - E bulll )/A (37) 


We now outline the procedure for doing the same calculation with ECT. Again we start with the 
same picture, an atom in the surface plane. We will only concern ourselves with up to next nearest 
neighbors in this example. Referring to equation (33), the perfect fee equivalent crystal has 12 nearest- 
neighbors and 6 next nearest-neighbors, therefore b^ = 12 and b 2 — 6. An atom in the free surface has 
lost four nearest-neighbors and one next nearest neighbor, thus we can rewrite equation (33) for the 
surface atom as 



An atom in the second plane keeps all of its nearest-neighbors and loses one next nearest-neighbor, 
therefore we have 


12 Rj (2) 


e 


-«Ri(a) 


+ 


6RP(2)e' aR2 ^e _R2 ^ /A 



(39) 


where now the arguments of R refer to the equivalent atom in each plane and noting that R 2 is simply 
a geometrical factor times R x and R’ is known since we are specifying the geometry of the defect as 
are a and A, which are given in the ECT papers. Thus, we are simply left with the problem of solving 
each transcendental equation for Rj for each nonequivalent atom in the surface region. Once this is 
done we simply substitute the values into equation (34) and evaluate the energy, where, for the rigid 
crystal, only the first term contributes. It is conceivable that for high degrees of symmetry and trun- 
cating at next-nearest neighbors that the calculation can be done on a programmable hand calculator. 

Let’s consider how many terms must be kept for the surface energy calculation. Note that the epu 
time consuming step in EAM is overlapping the densities, and in ECT it is solving the transcendental 
equation. Table I shows the contributions to the surface energy for EAM and ECT for a Ni(lOO) surface. 
In either case, three planes are sufficient to evaluate the energy. Table II shows the ECT values for an 
Al(210) surface. Since this is a high index plane, lower lying planes can have neighbors in the surface 
layer and thus more planes must be kept in the calculation. 

In table III we show a comparison between ECT, EAM, first-principles and experiment for the surface 
energies. ECT gives quantitatively accurate agreement with first-principles calculations (typically 
<10 percent disagreement) whereas EAM gives about 40 percent disagreement. In addition, ECT can be 
simply used for any crystal structure whereas EAM has been primarily applied to fee metals although it 
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has been used for bcc and hep in some cases. In table IV we show a comparison for surface relaxation 
and again ECT reproduces the magnitudes and trends quite well. This is a very stringent test, since the 
energy differences are very small in this case. In table V we show the results of using the more general 
form of ECT (eq. (36)) applied to Si. 

It is evident that EAM underestimates the surface energy. Our opinion is that these effects are 
caused by the fact that overlapping atomic densities does not include relaxation of the electron gas. At a 
free surface there are large deviations from equilibrium using simple overlap of free atom densities, where- 
as for an internal interface such as a grain boundary simple overlap is a less severe approximation and the 
energies obtained may be quite good. ECT has in some sense, electronic relaxation built into the proce- 
dure of finding the equivalent crystal and thus avoids this difficulty. Since the energy depends varia- 
tionally on errors in the electron density, errors in the density only cause second order errors in the 
energy. In this section, it is hoped that we have communicated that the application of EAM and ECT to 
problems of interest can be relatively simple. 


FRICTION AND OTHER STUFF 

Although these proceedings involve friction, we are going to avoid the real issue of how do you model 
the loss processes and instead present results of another problem which is a bound on the friction process, 
i.e., sliding one single crystal surface over another and also discuss a recent first-principles calculation of 
slip. 

We start with a general discussion of the adhesion by referring back to figure 2. Here we show a 
solid which is separated perdendicular to some plane. If we calculate the total energy at each separation, 
the adhesive binding energy is defined as 

E AD (a) = (e (a) - E(~))/2(2 A) ( 40 ) 

where a is the separation between surfaces and A is the cross-sectional area. Figure 5(a) shows what 
general shape one would expect the curves to have. Taking the derivative of this curve, the tensile stress 
is the curve shown in figure 5(b). Figure 6 shows the results of a first-principles, jellium calculation 
performed by Ferrante and Smith (ref. 29), solving equations (1) to (5) for this approximation. We can 
see that the curves for different metals in contact have large binding energies and have the same features 
as the same metal contacts. Figure 7 shows that all of these curves scale onto one binding energy rela- 
tion: the UBER described earlier. We now show adhesion calculations using ECT. The adhesion cal- 
culation for ECT or EAM is essentially as simple as the surface energy calculation. The only difference is 
that the other half space is present and its position is changed. Figure 8(a) shows (refs. 30 and 31) the 
results of an adhesion calculation for the (111) surface of a number of fee metals and figures 8(b) and (c) 
show the adhesive energy curves for the (110) surface of two bcc metals. The corresponding scaled 
binding energy curves are presented in figures 9(a) and (b) along with a fit to the Rydberg function. We 
can see that the curves scale. Figure 10 gives the force as a function of separation for (110) iron and 
tungsten. ECT can be formulated completely in terms of force, there-fore these calculations can be done 
directly. ECT calculations are trivial compared to the jellium or fully three-dimensional first-principles 
results presented. 

An interesting aspect of the adhesion calculation which was ignored in previous calculations relates to 
instabilities which can occur as surfaces approach on another called “avalanche.” Studies of this topic 
were prompted by a paper by Pethica and Sutton (ref. 32). The authors using ECT (ref. 33) and planar 
surfaces extended their results. The assumption that adhesion could take place with the half spaces 
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essentially rigid with possibly some small relaxation is clearly not correct. We found, as did Pethica and 
Sutton, that at a certain separation the surfaces will snap together depending on the thickness of the 
metals in contact. In retrospect this result is not too surprising since the spring constant for the restoring 
force depends on the length of the solid. Figure 11 shows the energy for moving one surface layer as a 
function of separation. At a separation of 2.6 A there is a barrier to moving one layer from each surface 
whereas at 1.9 A this barrier is removed. Figure 12 shows the energy as a function of separation for 
relaxing a number of surface layers where the circles represent the “rigid adhesion” curves from our 
previous calculations. As the number of planes relaxed (length of the solid) is increased, “avalanche” 
occurs at larger separations. Figure 13 shows how the separation at which avalanche occurs depends on 
the number of layers relaxed and the rigid separation. This example demonstrates that these calcula- 
tional methods enable examination of phenomena which are difficult to observe experimentally. 

We now address some ECT results that model friction more closely, that of the slip of one half- 
space over another (ref. 34). We have used ECT to calculate the energy and force as a function of tan- 
gential position for slip on an fee (100) surface for a number of metals. This calculation was performed 
under zero applied load conditions, i.e., as a force tangential to the surface was applied the two half- 
spaces were allowed to breath at the separation plane in the direction normal to the (100) surface. Once 
this calculation was finished the results were fitted to a Fourier series for the energy with a modulating 
term due to the normal displacement given by 


where 


°i T ) = X) A i I 1 + z i + ft z i1 e * iH i( x >y) 

(41) 

1=0 


z i " ( z " z t)/ 1 i 

(42) 

H o( x >y) = 1 - Hjky) 

(43) 

Hi( X ,y) = cos (27rx/a) cos(2jry/a) 

(44) 

H 2 (x,y) = cos(4ffx/a) + cos (4?ry/a) - 2 H^y) 

(45) 


where a is the lattice parameter and A-, 6j, and \ y are fitting constants from the ECT calculations and 
are given in reference 34. Equation (41), which can be thought of as an extension of the Frenkel model 
(ref. 1), is an excellent representation of the numerical results, with average deviation of order 10' 4 eV 
per surface atom for all metals treated. This is illustrated in figures 13 and 14 where the numerical 
results of equations (41) to (45) are plotted as solid curves. Figure 13 shows ideal adhesion or cleavage 
as a func-tion of z at fixed (x,y) for the Ag (100) interface. The bottom curve corresponds to (x,y) — 
(0,0), the middle to (a/4,0) and the top to (a/2,0). 

Equations (41) to (45) and figure 14 apply to arbitrary slip directions where the repeat distance can 
be many times the lattice constant a. Figure 14 shows zero-load results for two illustrative slip 
directions, the lower curve shows results for slip along the x-axis and in the upper an angle of arctan 
(3/4) relative to the x-axis. Figure 14 also shows the scaled results for all four metals examined where 
the energy was divided by the surface energy, A 1? and the distance by the lattice parameter, a. Thus 
there appears to be a universal form for ideal slip. This universality leads to useful rules of thumb for 
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slip, analogous to the Griffith criterion for cracks. Consequently, we Find that the average energy barrier 
to slip is E m = 0.36 Aj and the average maximum slip stress is < 7 m = 1.45 Aj/a. 

Figure 15 provides an overall picture of the stress profile for an Ag interface sliding on the (100) sur- 
face at zero load. The size of the arrow head and the length of the vector are both proportional to the 
magnitude of the slip stress. In table VI we show the numerical results for the height of the energy bar- 
rier and the maximum shear stress in the direction of minimum resistance to slip along with a comparison 
with whisker strengths. Although the theoretical results are larger, they are of the same order of 
magnitude as the experimental results. Since in the fee system slip occurs on (111) planes, the results are 
in reasonable agreement. The ability to map the entire energy surface and then analyze the results for 
other relationships shows the power of the semi-empirical methods. What we have presented here is not 
friction, however, since there are no energy dissipation mechanisms in the calculation. In addition, at 
this stage of the calculation we have ignored relaxation and any dynamic effects that will occur during 
sliding. However, the minimum stress to initiate sliding is a number of interest and the fact that scaling 
was found may mean that there are some underlying simplicities in the mechanisms. 

We now discuss a first-principles calculation similar to the above of sliding at a Pd-graphite inter- 
face by Zhong and Tomanek (ref. 35). They calculate the force to slide a mono-layer of Pd in registry 
with the graphite surface. Instead of a zero load, they calculate the sliding force as a function of load. 
Figure 16 shows their results for the position of the potential energy barrier and tangential force as a 
function of sliding distance in a given direction from hollow to bridge sites. The results are obtained from 
solving the Kohn-Sham equations (eqs. (1) to (5)) for the situation described. The effects of the load are 
obtained by calculating the adhesive energy, i.e., the energy as a function of distance in the normal 
direction. No specific mechanism is presented for energy dissipation and the assumption is made that all 
of the energy gained in being pushed up the hill is efficiently transferred into creating surface phonons 
and electron-hole pairs. Zhong and Tomanek then assume that the energy dissipated is given by 


AE f = (fj)Ax 


(46) 


where (f^) is the average friction force in sliding a distance Ax. The energy dissipated is estimated 
from the maximum barrier height in a given direction (fig. 16), and thus AEj = AV max anc ^ thus the 
average friction coefficient is given by the average friction force divided by the external load as 

. ■ <ffVf,x. - AV m „/(f ext Ax) («) 

A plot of this result is given in figure 17. The apparent anomaly for low load is a result of the 
minimum in the potential energy switching from a hollow to a bridge site under load. They then argue 
that the value obtained /z 0.02 is in agreement with experiment. 

One final topic we should address is ceramic and ionic compounds, since tribology has a wide range of 
materials involved in sliding contact and ceramic compounds are becoming increasingly important in 
newer high-temperature engines. At present pair-potentials (refs. 36 and 37) are the forms most often 
used to model interfaces and other defects. This approach suffers from the same difficulties previously 
described with the added complexity of having the possibility of having net charges at interfaces. 
Development of new approaches to treat these problems is definitely a goal for the future. 

First-principles calculations are presently being used to treat the problem of ceramic interfaces with 
their advantage of not being material dependent, but with the limitation that super-cells are used. As an 
example of such calculations we show a super-cell (fig. 19) structure used to calculate the interfacial 
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energy for a diamond C/sphalerate Boron-Nitride interface (ref. 38) using a linear-muffin-tin-orbital-approach 
(LMTO). The calculated adhesive energies range from 0.89 to 0.95 eV/unit cell area, depending on what 
structure is assumed for the interface. These calculations did not include any relaxation of the structure. 
These structural relaxations can be done if the changes are mainly local, which is probably the case. The 
first-principles calculations give promise for providing information concerning interfacial energies at least 
for the materials combinations of interest in tribology. 

At present there are no semi-empirical approaches for nonmetals. Recently, we have proposed (refs. 39 
and 40) an expression for the total energy for situations with charge transfer such as ionic solids with the 
hope that a method could be developed that would treat ceramic materials. The expression we propose is 
given by 


E(r) = -C (l + a*) e -a * - a£Z 2 f ( r)/R ( 48 ) 

where the first term allows for the possibility of both covalent bonding and the repulsion under bonding 
and the second the coulombic charge transfer, where C is the well depth of the Rydberg part, a is now 
(R - Rg)/1 and Rg is the distance at which the minimum occurs for the Rydberg part, £Z is the 
charge transfer, a is the modeling constant, and f(R) is a function that accounts for the crossing of the 
energy from two ions to two neutrals at some distance. This result can be contrasted with a Born-Mayer 
potential (ref. 1) which is a similar expression, but has a repulsive exponential core and thus does not 
allow for the possibility of covalent bonding. Figure 20 shows a fit of this function to the potential 
energy curves from first-principles calculations for three ionic molecules, A1C1, A1F, and LiF (a = 1, for 
this case). In addition to the fits being quite good the values obtained for the charge transfer, the 
relative covalent contribution, and the spectroscopic constants obtained from the global fits are quite 
reasonable. The hope in pursuing these studies is that a semi-empirical method can be developed similar 
to EAM or ECT in order to treat nonmetal and metal-nonmetal defects. 


CONCLUDING REMARKS 

In this presentation we have tried to give a description of the various computational techniques 
presently being used to approach problems in computational material science. We have also indicated 
that such calculational techniques now are sufficiently accurate and along with improvements in com- 
puter speed enable serious consideration for attacking problems in tribology for certain classes of mate- 
rials. At present it seems unlikely that the configurational energy minimization and dynamics needs for 
modeling tribological problems will be amenable to first-principles approaches because of the complexity 
of the problem and the computer time needed to treat each configuration. A more likely scenario will be 
the development of semi-empirical approaches which give accurate values for defect energies as compared 
to static configuration first-principles calculations and are then used to treat configurational energy 
minimization and dynamic effects. 
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TABLE I. — PLANAR CONTRIBUTION 
TO THE SURFACE ENERGY USING 
EAM AND ECT FOR Ni(100) 


Plane 

EAM 

ECT 

i 

0.290 eV/atom 

0.654 eV/atom 

2 

0.0182 

1.25X10 -6 


TABLE II. — PLANAR CONTRIBUTION 
TO THE SURFACE ENERGY USING 
ECT FOR A1 (210) 


Plane 

ECT 

i 

1.25 eV/atom 

2 

0.404 

3 

0.055 

4 

1.25x 10 6 
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TABLE III. — RIGID AND RELAXED SURFACE ENERGIES IN 
erg/cm 2 FROM ECT, FIRST-PRINCIPLES AND EAM 


(OTHER) CALCULATIONS IN REFERENCE 33 


Element 

Crystal face 

ECT rigid 

ECT relaxed 

LDA rigid 

Other 

Cu 


1830 

1780 

2100 

1170 


■i 

2380 

2320 

2300 

1280 



2270 

2210 


1400 

Ag 

(in) 

1270 

1230 


620 


(100) 

1630 

1600 

1650 

705 


(110) 

1540 

1510 


770 

Ni 

(111) 

2400 

2320 


1450 


(100) 

3120 

3040 

3050 

1580 


(110) 

2980 

2910 


1730 

A1 

(111) 

920 

860 




(100) 

1290 

1220 




(110) 

1310 

1230 

1100 


Fe 

(110) 

1820 





(100) 

3490 


3100 

1693 

W 

(110) 

3330 





(100) 

5880 


5200 

2926 
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TABLE IV. — PERCENTAGE CHANGES IN INTERLAYER 
SPACING DUE TO RELAXATION, REFERENCES FOR 


EACH VALUE ARE GIVEN IN REFERENCE 33 


Element 

Ad m - 

ECT, 

% 

EAM, 

% 

Experiment, 

% 

Technique 

Cu(110) 

Ad 12 

-7.7 

-4.9 

-8.5+0.6 

LEED 





-7.5+1.5 

Ion scattering 


Ad 33 

+3.4 

+0.2 

+2.3+0.8 

LEED 





+2. 5+1. 5 

Ion scattering 

Cu(100) 

Ad l3 

-3.7 

-1.4 

-2. 1+1.7 

LEED 





-1. 1+0.4 




*d 33 

+ 1.9 

-0.3 

+0.45+1.7 







+ 1. 7+0.6 

1 

r 

Cu(lll) 

Ad i3 

-3.1 

-1.40 

— 0.7+0.5 

LEED 


Ad as 

+ 1.9 

-0.05 




Ag(llO) 

Ad 13 

-6.0 

-5.7 

-5.7 

LEED 





— 7.8+2.S 

Ion scattering 


*d a3 

+2.8 

+0.3 

+2.2 

LEED 





+4.3+2. 5 

Ion scattering 

Ag(lOO) 

Ad ia 

-3.0 

-1.90 






+1.7 

-0.05 




Ag(Ul) 

Ad, a 

-2.6 

-1.30 





^d as 

+ 1.6 

-0.04 




Ni(110) 

Ad ia 

-7.6 

-4.87 

-8.7+0.5 

LEED 





— 9.0+1.0 

Ion scattering 


Ad a3 

+3.4 

-0.57 

+3.0+0.6 

LEED 





+3.5+1. 5 

Ion scattering 

Ni(100) 

Ad la 

-3.7 

-0.002 

— 3.2+0.5 

Ion scattering 


*d a3 

+2.0 

- .001 




Ni(Ul) 

Ad ia 

-3.1 

-0.05 

-1.2+1. 2 

LEED 


^d a3 

+ 1.9 

+ .00 




AI(llO) 

Ad la 

-10.4 

-10.4 

—8. 6+0. 8 







— 8.5+1.0 




Ad as 

+4.7 

+3.1 

+5.0+1. 1 







+5. 5+1.1 

1 


Al(lOO) 

Ad 12 

-4.9 






^d a3 

+ 1.8 





AI(lll) 

Ad 12 

-3.9 


+0.9+0.7 

LEED 


^d a3 

+2.5 
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TABLE V. — SURFACE ENERGIES IN ergs/cm 2 FOR Si(100) 


[The (2x1 


values are for symmetric dimer formation. 


Surface 

ECT 

Reference 

Reference 

Reference 

Reference 



a 

b 

c 

d 

Ideal 


mm 




lxl 

1 

IB 




2x1 

1550 


1590 

1610 



a. M.T. Yin and M. Cohen, Phys. Rev. B 24, 2303 (1981). 

b. K.C. Pandey in Proceedings of the seventh International 

Conference on Physics of Semiconductors, eds. D.F. Chadi 
and W.A. Harrison (Springer-Verlag, New York, 1985), p. 55. 

c. N. Roberts and R.J. Needs, J. Phys. Cond. Mat. 1, 3139 

(1989). 

d. M.I. Baskes, J.S. Nelson, and A.F. Wright, Phys. Rev. B 40, 

6085 (1989). 


TABLE. VI — ENERGY BARRIER TO SLIP 
E m (eV/atom), MAXIMUM SLIP STRESS 
a m (10 10 dynes/cm 2 ) AND MEASURED 
[a] MAXIMUM WHISKER STRENGTH 
o (IQ 10 dynes/cm 2 ) 


Element 

mm 

■hvstu 

a 

Ag 

0.292 

6.19 

3.32 

Ni 

.411 

12.3 

3.94 

Cu 

.334 

9.36 

3.64 

A1 

.225 

4.34 

— 


a. Whiskers, ed. J. Gordon Cook (Mills and 
Boon Limited, London, 1972) pp. 46-55. 
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Electron Density (e/nrn^) 



Figure 1.— A 4° rigid twist angle grain boundary for an fee (100) 
surface. 



Figure 3.— A grain boundary represented by a dip In the jellium 
density where the square lines are the jellium density and the 
Irregular the relaxed electron density 



POSITION y (nm) 

Figure 2. — A bimetallic AI-Mg Interface where the square lines 
are the Jellium density and the irregular lines the relaxed the 
electron density. 



0 2.0 4.0 60 7.0 

SCALED SEPARATION a* 

Figure 4. — Scaled binding energy as a function of scaled 
separation for the four systems shown. 
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ADHESIVE ENERGY (erg/cm2; 



1 

(a) Binding energy versus separation. 



ADHESIVE BINDING ENERGIES 






SEPARATION a(nm) 


Figure 6.— Adhesive binding energy versus the separation, a, 
between the surfaces Indicated. 



Figure 7. Scaled adhesive binding energies versus scaled 
separation for the Interfaces indicated. 




0 

-.2 





0 1 2 5 N 


i , ANGSTROM 

(a) For the (111) surface of Al p Ni, and Ag. 



DISTANCE BETWEEN SURFACES, A 
(b) For the W(1 1 0) surface. 



DISTANCE BETWEEN SURFACES, A 


(c) For the Fe{1 1 0) surface. 

Figure 8.— Rigid adhesive binding energy versus 
separation using ECT. 


23 




U01V/NU 



(a) fcc(1 11) Al, Nf, Cu and Ag. (b) bcc(1 1 0) W and Fe. 


Figure 9. — Scaled adhesive energies from Figure 8. 



Figure 10. — Rigid adhesive binding force versus separation for Figure 1 1 . — Total energies as a function of symmetric movement 

W(1 1 0) and Fe(1 1 0) interfaces. of only the surface atomic layers for N 1(1 00) for two rigid 

separations. 
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Energy (eV /rurfac 


0 




Interfad&l Separation (A) 

Figure 14. Total cleavage energy as a function of interfacial 
separation for Ag(1 00) for different positions on the surface. 
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0 1 2 3 4 5 6 

RIGID INTERFACIAL SEPARATION d R (A) 


Figure 13.— Relaxed interfaclal separation d as a function of 
rigid (unrelaxed) separation, d R for different number of 
layers allowed to relax. 



Scaled Distance 


Figure 1 5.— Scaled total energy as a function of scaled slip 
distance for two different slip directions on (1 00) Interface 
of Nl, Cu, Ag, and Al. 





Figure 16.— A map of slip stresses for Ag crystals sliding on (100) 
surfaces at zero load. 



0,0 2.5 5.0 


x(A) 

(b) Friction force as a function of position of the same 
interface [35], 


Figure 17. — Interfacial potential energy and friction force. 







Friction Coefficient /x 



Figure 1 8.— Microscopic friction coefficient p as a function of the 
external force per atom f ext . 


DIAMOND C SPHALERITE BN 



[ 110 ] 


• c O N © B 

Figure 1 9. — Crystal structure of the C2(1 1 0) and BN(001 ) 1 xl 
superlattice [39], 
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